In a number of situations in a variety of fields, ranging from reconstruction of MR images, other medical image reconstructions, image deblurring, to temperature estimation in transient problems (including automotive applications) and vibrations and geophysics (such as the seismic inverse problem), the problem reduces to reconstructing an object from observations of interaction of the object with a physical system.
One illustrative but very practicality important, example of the solution of this type of problem is the reconstruction of MR images. Magnetic resonance (MR) imaging has great importance, both for clinical and for research applications, due to its safety and exibility. The MR imaging process, however, imposes a fundamental tradeoff between image quality and scan time. It is very important to reduce scan time, for a number of reasons. The primary issue is that MR is very sensitive to motion artifacts, so reducing scan times decreases the odds of obtaining a poor-quality image. Some parts of the body, notably the heart and lungs, undergo periodic motion; as a result, only sufficiently fast scans can be used to image these organs, or any other adjacent tissues. In addition, there are obvious issues of patient comfort as well as cost that arise from lengthy scans. Finally, if the basic image acquisition process were faster, even for a xed scan duration a higher signal to noise ratio could be obtained. This is important for many important modalities, such as fMRI, perfusion, diffusion or time-resolved angiography, where the physical process being imaged takes place over short time periods.
There is a need for better algorithms to result in faster (and hence better) MR. The image formation process in MR is quite unlike a conventional camera, and this in turn leads to a number of interesting computational challenges. In MR, an image is acquired as its Fourier transform (usually called “k-space”), and this acquisition takes place serially. Typically one row in k-space is acquired at a time, and each row takes tens of milliseconds to acquire (the details depend on the particular study being performed).
A particularly important technique for accelerating MR scans is called parallel imaging, which uses multiple receiver coils. According to a recent survey article “Parallel imaging is one of the most promising recent advances in MRI technology and has revolutionized MR imaging”. While parallel imaging requires solving a linear inverse problem, existing methods either do not incorporate spatial priors, or assume that the image is globally smooth. In computer vision, of course, Markov Random Fields (MRF's) are commonly used to encode better priors.
The parallel imaging problem is illustrated in FIG. 1a. FIG. 1a is a schematic of the pixel-wise aliasing process, for a single pair of aliasing pixels p and p′, for 2× acceleration using 3 coils. The aliased observations Y are obtained by a weighted sum of the aliasing pixels, weighted by coil sensitivity values S1. To simplify the figure, aliasing is shown in the horizontal direction. The scan time can be reduced in half by dropping every other column in k-space. The resulting image, when reconstructed from a conventional RF receiver coil, is an aliased image like those shown in the top row of FIG. 1a. Such an image is half the width of the original (unobserved) image shown in the bottom row. It is formed by multiplying the original image by a slowly varying function S(p) and then adding the left half of the resulting image to the right half. The multiplication by S(p) comes from the spatial response of the coil, while the addition of the two image halves is a consequence of dropping alternate k-space columns. Thus, each pixel pin the aliased image is the weighted sum of two pixels (p,p′) in the original image, where p′ is in the same row but half the width of the image away. The weights come from the function S.
If only had a single aliased image were available, it would be impossible to reconstruct the original image. However, we multiple coils can be used, each of which has a different S, without increasing the scan time. In the above example, there is a simple linear system for each pair of pixels p, p′:
                              [                                                                                          Y                    1                                    ⁡                                      (                    p                    )                                                                                                                                            Y                    2                                    ⁡                                      (                    p                    )                                                                                                                                            Y                    3                                    ⁡                                      (                    p                    )                                                                                ]                =                                            [                                                                                                                  S                        1                                            ⁡                                              (                        p                        )                                                                                                                                                S                        1                                            ⁡                                              (                                                  p                          ′                                                )                                                                                                                                                                                S                        2                                            ⁡                                              (                        p                        )                                                                                                                                                S                        2                                            ⁡                                              (                                                  p                          ′                                                )                                                                                                                                                                                S                        3                                            ⁡                                              (                        p                        )                                                                                                                                                S                        3                                            ⁡                                              (                                                  p                          ′                                                )                                                                                                        ]                        ⁡                          [                                                                                          X                      ⁡                                              (                        p                        )                                                                                                                                                        X                      ⁡                                              (                                                  p                          ′                                                )                                                                                                        ]                                .                                    (        1        )            
The functions S are called sensitivity maps, and can be computed in advance. In addition, the sum over the coils is usually normalized:
            ∑              l        =        1            L        ⁢                  ⁢                  S        i        2            ⁡              (        p        )              ≈  1.
Formally, assuming Cartesian sampling in k-space, an (unobserved) true image X is desired, which must be reconstructed from the sensitivity maps S1, . . . S1, . . . SL and the coil outputs Y1, . . . Y1, . . . YL. While X and S1 are of size NxM, Y1 is of size
            N      R        ⁢          M      x        ,where the acceleration factor is R. By stacking the receiver coil outputs into the column vector y and the image into x, the image formation process can be expressed asy=Ex  (2)                where the encoding matrix E contains the sensitivity maps.        
Least squares can be used to compute the maximum likelihood estimate from equation (2){circumflex over (x)}SENSE=(EHE)−1EHy  (3)
This is an MR reconstruction algorithm called SENSE, which has been widely applied in research and clinical systems. Tikhonov-style regularization methods are also commonly used to handle the ill-conditioning of H. The most general form of these methods is
                                          x            ^                    regSENSE                =                                            arg              ⁢                                                                                ⁢                                                                              ⁢              min                        x                    ⁢                      {                                                                                                  Ex                    -                    y                                                                    2                            +                                                μ                  2                                ⁢                                                                                                A                      ⁡                                              (                                                  x                          -                                                      x                            T                                                                          )                                                                                                  2                                                      }                                              (        4        )            
Here, the second term penalizes that are far from a prior reference image xr, and that are not spatially smooth according to the regularizer matrix A. The regularization parameter μ controls the tradeoff between aliasing and local noise in the output. Note that since A is a linear operator, global spatial smoothness is imposed. There is a simple closed-form solution for this equation.
FIG. 1a is typical of the parallel imaging problem, for 2× acceleration with 3 coils. The unobserved image is shown at bottom, and the aliased image observed from each coil is shown at top. Each pixel in an aliased image is the weighted sum of two pixels in the unobserved image. The weights are different for each pixel and each aliased image. The linear system that relates the unobserved image to the observations is given in equation (1).
Every clinical system, and almost every research paper in MR, computes some variant of x^reg. The original paper on SENSE did not perform regularization (i.e., μ=0). For numerical reasons, authors proposed variants of equation (4), typically using the identity matrix for A. While this does not impose spatial smoothness, it has significant computational advantages in terms of the structure of H. Even the most general variants of this method, however, can only impose spatial smoothness globally, and give poor performance for acceleration factors significantly greater than 2.
Other reconstructions of an object from observations of interaction of the object with a physical system have similar detailed descriptions.